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In the present paper we consider an optical system with a x^-type nonlinearity and unspecified 
PT-symmetric potential functions. Considering this as an inverse problem and positing a family 
of exact solutions in terms of cnoidal functions, we solve for the resulting potential functions in a 
way that ensures the potentials obey the requirements of PT-symmetry. We then focus on case 
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I. INTRODUCTION 

Recently a great deal of attention has been devoted to the investigation of quantum and classical systems with PT- 
symmetric non-Hermitian Hamiltonians. It was shown in the seminal paper [1] that such type of Hamiltonians can 
have real eigenvalues. Due to the analogy between the Schrodinger equation and the paraxial wave equation in optics, 
this result has applications in optics; this is one among numerous other areas studied over the past decade. For optical 
beams PT-symmetry imposes the condition on the complex refractive index n(x) = n r (x) + irii(x): even in space 
for the real part of the refractive index n r (x) = n r (—x) and odd in space for the imaginary part rii(x) = —rii(—x). 
Recently effects of PT-symmetry have been observed in optical experiments [2]. However, optics is certainly not the 
sole area where PT-symmetric applications have recently emerged. More specifically, a mechanical system realizing 
PT-symmetry has been proposed and realized in [3j, while a major thrust of efforts has focused on the context of 
electronic circuits; see e.g. the original realization of |4] and the more recent review of this activity in [5]. Additionally, 
we note that further intriguing realizations of PT-symmetry have also emerged e.g. in the realm of whispering-gallery 
microcavities [Bj. While many of these experimental realizations have been chiefly explored at the level of linear 
dynamics, the intrinsic nonlinearity of optical systems [2j and the potential nonlinearity also of electrical ones (e.g. 
in the form of of a PT-symmetric dimer of Van-der-Pol oscillators Q) have prompted a considerable amount of work 
at the interface of nonlinearity and PT-symmetry. 

Nonlinearity leads to new effects in the PT-symmetric systems, such as solitons (and vortices) in continuous [8j[9 
and discrete nonlinear optical media with PT-symmetric potentials m, gap solitons in media with PT-symmetric 
periodic potentials m , non-reciprocity, instabilities and nonlinear PT-transitions in PT-symmetric (nonlinear) cou¬ 
plers mm, as well as the smoothing of the spectral singularity in transmission [18], among many others. Important 
applications to nonlinear plasmonic systems and metamaterials are under recent investigation H3H20. Additional 
developments are connected with nonlinear PT-symmetric lattices [2TH24] and PT-symmetry management [25l [27 
etc. 

In nonlinear wave equations with PT-symmetric terms, solitonic solutions can exist as was shown recently e.g. in 
the works of [51 1281131] . In the case of an NLS system with inhomogeneous in space loss/gain parameters such solutions 
were found for linear PT-potentials in the work [H EH E3] , and for nonlinear PT-potentials in the works [28l 129] . 
for a nonlocal NLS equation in [34] and for a cubic- quintic model in [35] . Naturally, it is also of interest to find 
exact solutions for solitons in other physically important systems. Recent numerical simulations of the x^ system 
with PT-symmetric potential showed the existence of stable bright solitons [36], as well as of gap solitons in periodic 
PT-symmetric potentials EH- Discrete PT-symmetric systems with quadratic nonlinearity were also explored at the 
level of oligomer systems [38] . 

In the present work we will study the x^ system with PT-symmetric potentials, describing wave processes in 
quadratically nonlinear media with spatially distributed gain/loss parameters. Such systems are, in principle, of 
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interest for nonlinear optics and potentially even for atomic-molecular Bose-Einstein condensates trapped in complex 
potentials (see, for example, a realization of imaginary potential in [39] )• In this paper, we find exact solitary and 
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periodic solutions of the x ^ system. We begin in Section |n| by outlining the mathematical model. In Section 
we derive exact solutions in terms of the Jacobi elliptic cnoidal function, and we present special case solutions for 
which we later examine stability properties. We essentially follow an inverse-function approach, somewhat reminiscent 
of [00] to obtain such exact solutions for suitably tailored PT-symmetric potentials in the presence of the quadratic 
nonlinearity. In Section |IV| we present the stability analysis of the obtained solutions as a function of the solution 
parameters (such as their amplitudes), and we show the results of propagation of the solutions in the (analogous to 
“time”) variable z. Finally, in Section [V] we make our concluding statements and present a number of possibilities for 
future work. 


II. THE MODEL 

Let us consider the y 2 system describing the first harmonic (FH) and second harmonic (SH) propagation in quadrat- 
ically nonlinear media with PP-symmetric potentials as follows: 

iu z + di u xx + Vi(x)u + iW\(x)u = u*v (1) 

iv z + d 2 v xx + kv + V 2 (x)v + iW 2 (x)v = u 2 

where V\^ 2 {x) are even functions of x corresponding to real parts of the refraction index, and Wi, 2 (x) are odd functions 
of x pertaining to imaginary parts thereof. W \,2 describe the inhomogeneous in space gain/loss. Seeking standing 
waves in the form: u(x, z) = U{x)e~ luJZ and v(x, z) = V(x)e~ 2luJZ we obtain the system 

ujU + dxU^ + Vi^U + iWi(x)U = U*V ( 2 ) 

aV + d 2 V xx + V 2 (x)V + iW 2 (x)V = U 2 

where a = 2cj-\-k. It is useful to introduce the amplitude-phase decomposition in the form U(x) = pi{x)e ie ^ x \ V(pc) = 
p 2 (x)e 2l0 ( x \ Solitonic solutions for V = W = 0 have been reported e.g. in [4lti45] . cnoidal wave solutions in [46] . 
and solitons in the conservative 2D system with a potential V ^ 0, W = 0 were considered recently in [47] . For 

a review of solitary wave dynamics in quadratic systems see e.g. [48] . 

Assuming that p 2 ,0 are real, and that pi is either real or purely imaginary, we obtain the system 

P*iP 2 = upx + dxp ltXX - dipi{0 x ) 2 + V 1 (x)p 1 (3) 

Pi = a P 2 + d,2P2,xx ~ 4^2/32(#;r) 2 + ^(x)p 2 
Wj(x)p 2 j = - jdj{p 2 0 x ) x 


for j e {1,2}. 


III. SOLUTIONS IN TERMS OF THE CNOIDAL FUNCTION 

We begin by writing pi, 2 = Fi, 2 (y) and 0 X = G(y) for y = cn(rx, k). This gives the system 

Ft{y)F 2 {y) = uF^y) + d^T^y) - d 1 Fi{y)G 2 {y) + V^F^y) (4) 

F 2 {y) = aF 2 (y) + d 2 r 2 T 2 (y) - M 2 F 2 (y)G 2 (y) + V 2 (x)F 2 (y) 
w j( x )Fj(y) = jrd j da(rx,k)sn(rx,k)(2Fj(y)G(y) + G'(y)F j (y)) 

for 

I'M = y(2k 2 - 1 - 2 k 2 y 2 )F’M + (1 - y 2 ){k 2 y 2 + 1 - k 2 )F'p y ) (5) 

with j G {1,2} and where the primes denote differentiation with respect to y. Notice that by writing Q in terms 
of both x and y = cn(rx, k) we avoid restrictions on the domain which would be applicable if we composed with an 
inverse function. 

To find exact solutions, we apply the reverse engineering approach ^30^5 ^is type of technique was applied much 
earlier in order to obtain exact traveling wave solutions in dynamical lattices m • Our general strategy is to first 
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specify the form of the functions Fi^(y),G(y) and then use © to solve for appropriate potentials Vi^{x) : Wi^{x). 
Thus, we rewrite the above equations as: 


Vi(x) 

V 2 (x) 


Wj(x) 


Fi(y) a „2 r i(y) , . ^2 


Fl (y) F2{y) - dir F l{ y) 
F 2 (y)-d 2 r 2 T 2 (y) 


+ d\G ( y ) — LO 


F 2 (y) 


jrdjdn(rx,k)sn(rx,k) I 2 


4 d 2 G 2 (y) - cf 
Fj(y)G(y) 


Fj(y) 


G'{y) 


(6) 


for y = cn (rx,k) and j G {1,2}. In each of the following subsections we make specific choices of Fi ?2 and G in 
such a way that the resulting potentials Vi j2 , VFi j2 in © do not contain terms with denominators (that may lead to 
singularities) and they also obey the requirements of PT-symmetry (Vq 2 even functions of x and W\^ odd functions 
of x). In other words, we require at least that the conditions 


Fi(y) 

1 r a (y) 

(7) 

F 2 (y ) 

| (F 2 (y) — d 2 r 2 T 2 (y)) 

(8) 

Fj(y) 

1 7(y)G(y) 

(9) 

for j = 1 , 2 are met for any choices of Fi j2 , G that 

we specify. 



A. Polynomial Functions 

Consider Fi^(y) in the form of generalized polynomials in y — cn(ra, k ) 


k 2 


Fi(y) = i {0 ’ 1] £ C n y n , F 2 (y) = D " 


( 10 ) 


with coefficients C n , D m G IR, integer indexing bounds si, &i, > 0 with k\ > si, & 2 > s 2 , and C Sl , Ck 1 , D S2 ,Dk 2 7 ^ 

0. Notice that we restrict our attention here to polynomials with at least two terms. The case of a monomial type 
solution will be included in the next section where we consider a more general class of power functions. Recall that, 
in the derivation of ©> pi is required to be either real or purely imaginary. To show this in ( [To| ) we have included 
an optional multiple of i in the definition of F\. In the following subsections we outline the process of solving for 
Vi j2 , W h2 - We separate into two cases which are convenient based on the resulting maximal power of the polynomial 
conditions ©-©• 


1. Cnoidal parameter k ^ 0 

To proceed in solving for V\(x) using © and assuming ( [To] ) we must satisfy condition ©. That is, we must have 

that Fi(y) is a factor of the polynomial Ti (y). For k 7 ^ 0, Iiu/) will have maximal power k± + 2 so that 0 amounts 

to the condition 

fci+2 k\ 

T, ((n + 2 ){n + 1)(1 - k 2 )C n+2 + n 2 ( 2 k 2 - 1 )C n - [n - 2)(n - 1 )k 2 C n - 2 ) y n = (any 2 + fay + 71) ^ C n y n 

n=si —2 n=s\ 

(n) 

for some aq, /?i, 71 G IR with ol\ 7 ^ 0. For convenience we use the convention that Gj = 0 for any index j 0 {si,..., /q}. 
Equating the coefficients of then gives 

(n 2 ( 2 k 2 - 1 ) - 71) C n = (an + (n - 2)(n - l)/c 2 ) C„_ 2 + AC„_i - (n + 2)(n + 1)(1 - fc 2 )C n+2 (12) 

«i = -An (An + l)k 2 , si(si - 1)(1 - k 2 ) = 0 , (si + l)si(l - k 2 )C Sl+ i = 0 

where the first equation is a recursion relation that holds for n € {si,..., An + 1 } and the latter three equations are 
obtained from equating the coefficients of the y kl+2 , y Sl ~ 2 , y Sl_1 terms in < 00 , respectively. Note that the latter 
equations have made use of the conditions C Sl 7 ^ 0 and Cj = 0 for j 0 {si,..., ki}. 
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From the first equation of the latter three in (12), we now have that the coefficient oq is determined by the highest 
degree chosen for F\. The latter two equations in ( 12 ) then give us a starting point for finding more specific solutions. 
That is, we can look for solutions with k = 1 in terms of cn(rx, 1) = sech(rx) and these latter two equations are 
satisfied. Alternatively, we can look for solutions with k ^ 1 in which case the latter equations of ( 12 ) give that either 


<§1 =0 so that the first term in the polynomial F\ is required to be a constant, or the first term is required to be 
of degree s i = 1 with coefficient C Sl +% = G 2 = 0 . We will proceed here to outline the general solution for cnoidal 
parameter k 7 ^ 0. Later towards the end of this section we will focus primarily on k = 1 for the special case where 
Fi ? 2 are quadratic polynomials. 

The recursive equations in ( 12 ) give us Aq — si + 2 conditions for the Aq — si+3 constants {/3i, 71 , C Sl , G Sl +i,..., G/^}. 
Later we will choose f3\ and then use conditions to solve for the coefficients of F\ and also 71 . Now that 0 is 
satisfied by imposing ( 12 ), we have V\(x) via 0 as 


Vi(x) = ±F 2 (y) - d 1 r 2 {a 1 y 2 + foy + 71 ) + d 1 G 2 {y) - cj 


(13) 


with y = cn(rx, k) as usual. The plus sign in ( [T3| ) applies to F\ real (using i° = 1 in ( |Tq| ) ) and the minus sign applies 
to Fi purely complex (using i 1 = i in (|Tq|)) . Since the cnoidal function is an even function of x, V\(x) is an even 
function of x so that this potentia l is compatible with PT-symmetry. Notice that V\(x) may not be a polynomial if 
G 2 is not a polynomial. In Section [LV| we will specify choices for the G function and we will choose c u so that V\ —» 0 
as x —>> 00 . 

To solve for V 2 {pc) we must have that F 2 (y) statisfies condition One way to proceed is to require that 


F?(y) = F 2 (y)P(y) 


(14) 


where P(y) is a polynomial in y. Then, similar to the F\ case, we may also impose that F 2 (y) | T 2 (y) so that 
T 2 (y) = (ol 2 V 2 + fcy + 72)^2 (y) for some constants « 2 , ^ 2,72 G R with a 2 7 ^ 0 . Using a similar procedure as in the Fi 
case, now equations ( [l 2 | must hold after performing the replacements n m, C D and in the subscripts 1 —>> 2. 
Using this inversion of equation (12), now the coefficient <a 2 is determined by the highest degree of the polynomial 
F 2 . Since k 7 ^ 1 here the inversion of the latter two equations in either requires us to take s 2 = 0 so that the 
first term in the polynomial F 2 is required to be a constant, or alternately the first term is required to be of degree 
S 2 = 1 with coefficient C S2 +1 = C 2 = 0 . Also, the — s 2 + 3 constants {/3 2 ,7 2 , D S2 , D S 2 +i 5 ...,^ 2 } are required to 
satisfy the same recursive fc 2 — s 2 + 2 equations in (fl 2 l but with appropriate inversion described above. 

The most obvious choice in order to satisfy both (|14|) and the F 2 -version of ( 12 ) is to take F 2 (t/), Fi(y) as scalar 
multiples of each other. In other words, 


F 2 (y) = i^AF 1 (y) and P{y) = TiA 


(15) 


for some A E R/o- Since p 2 , U 2 are required to be real-valued the multiple of i in front is included only if it’s included 
in the definition of F\ in (10). Now we have the real-valued potential function 


V 2 (x) = - d 2 r 2 (a 2 y 2 + /3 2 y + 72 ) + 4 d 2 G 2 {y) - a 


(16) 


where later in specific examples we will choose a to be such that U 2 0 as x —* oo. 

Next we want to determine an appropriate form for the function 6 X = G(y) E R with y = cn (rx, k) that will satisfy 
0. Since F\,F 2 have been chosen to be scalar multiples of each other, if ([9| holds for j — 1 then it holds for j = 2 . 
So, we take 


G(y)=T(y)F 1 (y) 


(17) 


for a function T(y) and this gives via 0 


Wj(x ) = jrdjdn(rx , k)sn(rx, k) ( 2 F[(y)T{y) + G'(y)) 


(18) 


for y = cn(rx, k) and j E {1,2}. Since sn {rx,k) is an odd function of x and dn(rx, k) is even, VFi, VF 2 are odd 
functions of x as required by PT-symmetry as long as the quantity 2F 1 '(cn(rx, /c))T(cn(rx, k)) + G'(cn(rx, k)) is an 
even function of x. This is reasonable since cn(rx, k) is an even function of x. 

Now we have a complete solution of © given by Fi^ in ( p~ 0 | ), Ui ?2 in ( fl~3| ) and ( [Tg] ) , and VFi j2 in ( [18] ) , all under the 
conditions seen in ( [T 2 ] ) , ( fl5| ), To be more explicit, let us focus on details in the case where F\ is a quadratic 

function and k = 1 so that y = sech(ra). Consider of the form = Go + Giy + G 2 ^/ 2 for Go, G 2 ^ 0 . Then si = 0 
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and ki = 2 so that the latter two equations in ( 12 ) are satisfied. Now the remaining equations in ( 12 ) give aq = — 6 , 
71 = 0 and the conditions 


Ci =/?!(?0, C 2 = 


—6Cq + PiC\ 


c 3 = 0 = 


(<ai + 2)Ci + (3\ C 2 


(19) 


that we may solve for the four constants /?i, Co, Ci, C 2 . 

In the case that /3i = 0, (19) gives Ci = 0 and C 2 = —3Cq/2 so that combining with (15) we have 


Fi(y) = C 0 (i - g2/ 2 ), 


^ 2 ( 2 /) = AFi(y). 


Proceeding with G{y) as in © for any function T(y) we have by ([l3|, ( p~6| ) and ( p~8] ) that 

UOe) = ^4Co(l — ^2/ 2 ) + Qdir 2 y 2 + d 1 G 2 (y) — w 

r 2 (a;) = ^}(1 - ^y 2 ) + Qd 2 r 2 y 2 + 4d 2 G 2 (y) - a 
Wj{pc ) = jrdjtanh(rx)sech(rx) (— 6CoyT(y ) + G'(y)). 


( 20 ) 


( 21 ) 


Equations (20)-(21) give us a solution that we will refer to as the dark-dark soliton case. In Section IV we show the 
dark soliton shape, analyze the stability of the dark-dark soliton, and show plots over the propagation variable z for 
a specific choice of the function G and other parameters. 

We also consider here the case of f3\ ^ 0, for which (19) gives two possibilities for the coefficients of the polynomial 
F\. Then combined with (15) we have 


F 1 (y) = C 0 (l±V22y + 4y 2 ), F 2 (y) = AF 1 (y). 

Letting G(y) be as in © for some T(y) function we then obtain 

Vi(x) = ACo(l ± V22y + Ay 2 ) — dir 2 (±V22y — 6y 2 ) + diG 2 (y) — uj 
V 2 (x) = ^-(l±V22y + ^y 2 )-d 2 r 2 {±V22y-6y 2 )+M 2 G 2 (y)-a 
Wj[x ) = jrdjtanh(rx)sech(rx) ^2Co(8y ± V22)T(y) + G'(y)^ 


( 22 ) 

(23) 


for j G {1,2}. These solutions are quite interesting in their own right, as the one with the + sign corresponds to 
an antidark-antidark soliton setting of a pair of bright solitary waves on top of a non-vanishing background. On the 
other hand, the solution with the — sign is especially structurally complex, resembling a conglomeration of multiple 
-more specifically of 4- dark solitons. 


2. Cnoidal parameter k = 0 and y — cos(ra:) 


For Fi ?2 of the polynomial form (10) now consider the case of k = 0, or y = cos [rx). In solving for Vi(x) 
the polynomial condition analogous to (11) has maximal power k\. This roughly makes sense because when we 
differentiate a cos (rx) or sin(rx) the result is a function of the same overall power (in contrast to derivatives of 
sech(rx) and tanh(re), for example). The analogue of ( pTj ) in this case is 


k\ k\ 

V ((n + 2 )(n + l)C n+2 - n 2 C n ) y n = 71 V C nV r 

n=si — 2 n=s 1 


(24) 


for 71 7 ^ 0. Equating coefficients we get 


C„ = (" + 2 )(”+ 1 )0.w for nefo,. 


n z + 71 


7i = 


si(si - 1 ) = 0 , 


(si + l)siC Sl +i — 0 


(25) 


where the latter three equations came from equating the coefficients of the y kl , y Si 2 ,y Sl 1 terms in ( |24| ), respectively. 
As before, the first equation of the latter three in (25) shows that the coefficient 71 is determined in terms of the 
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maximal power k\ of the polynomial F\. The latter two equations in ( [25] ) then show that either s\ = 0 so Fi must 
have a constant term, or alternatively the first term is required to be of degree s i = 1 with coefficient C Sl +i = C 2 = 0 . 
The remaining recursive equations in (25j then give us k\ — s\ conditions for k± — s\ + 1 unknowns {C Sl ,..., C/^}. 
Choosing one of these coefficients will lead us to find the others. In solving for V 2 we have similar conditions to (25) 
for the constants 72 7 ^ 0 and {D S2 ,..., D^ 2 } where in (25) one should replace n —)> ra, C D and in the subscripts 
1—^2. We proceed in a similar way as in the k 7 ^ 0 case above, assuming the forms of F 2 ,G as seen in (14), (15), 
(p~7|) and finally we have 


Vi(x) = ±F 2 (y) - dir 2 7 i + diG 2 (y ) - u (26) 

V 2 (x) = ™L-d 2 r 2 l2 +M 2 G 2 (y)-a 

Wj(x) = jrdj sin(rx) (2F((y)T(y) + G'(y)) 


for j £ { 1 , 2 }. 

Let us focus on the details of the quadratic case where Fi = Cq + C\y + C 2 y 2 for Co, C 2 7 ^ 0 and sq = 0 , &q = 2. 

Then, we have 


(25) then gives 71 = —4, Co = — C 2 / 2 , and C\ = 0 . 

F x = C 0 (l - 2y 2 ), 


F 2 (y) = AFi(y). 


(27) 


We also have 72 = —4 by the V 2 analogue of ( |25| ) (see description above). Proceeding with G(y) as in for some 
function T(y) we have 


Vi(x) = ACo(l - 2y 2 ) + 4dir 2 + diG 2 (y) - cj 
V 2 {x) = ^(l-2y 2 )+4d 2 r 2 +4d 2 G 2 (y)-a 
Wj(x) = jrd j sm(rx)(-8C 0 yT(y) + G'(y)) 


(28) 


for j £ {1,2}. Equations (27 )-(|28[) give us a solution that we call the quadratic oscillatory case. In Section [TV] we 
show the shape of the solution, analyze the stability, and explore its dynamics over the evolution variable (z). 


B. Power Functions 


Next we take Fi(y) and F 2 (y) to be power functions 

F 1 (y)=iWCyP\ 


F 2 (y ) = DyP* 


(29) 


for pi,P 2 > 0 and C,D ^ 0. This special case considerably simplifies the relevant compatibility conditions. In 
particular, substituting (29) into <© and examining conditions 0-© we require that either k = 1 and 2 pi > p 2 , or 
Pi =P2 = 1. 


1. Cnoidal parameter k — 1 and y = sech(ra:) 


In the case of k = 1 we can have non-integer pi and p 2 ; this is in contrast to the polynomial case. For this case, we 
apply © and find 


Vi(x) = ±Dy P2 -d 1 r 2 p 1 (l-2y 2 )-d 1 r 2 p 1 ('p 1 -l){l-y 2 )+d 1 G 2 (y)-u 
V 2 (x) = ±^-y 2pi ~ P2 - d 2 r 2 p 2 ( 1 - 2 y 2 ) - d 2 r 2 p 2 (p 2 - 1)(1 - y 2 ) + 4d 2 G 2 (y) - a 


(30) 


Wj(x) = jrdj tanh(ra;) ( 2 pjG(y) + yG'{y)) 

for j G { 1 , 2 }. As long as ( 2 pjG(y) + yG'(y )) are even functions in x, which is easy to choose since y = sech(rx) is 
an even function of x, then W\ and W 2 are odd and compatible with the VT symmetry criterion. We refer to the 
solutions given in (29) and pQ|) as the bright-bright soliton case. In section IV we show the wave’s shape, analyze its 


stability and explore its direct numerical evolution. Note that the case pi = p 2 = 2 corresponds to solitonic solutions 
found by Karamzin-Sukhorukov in m, and 2pi = p 2 = 2 to solitonic solutions found by Menyuk et al. [44] . 
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2. Powers pi = p 2 = 1 


In the case pi = P 2 = 1, we have 

Vi(x) = ±Dy — d\r 2 (2k 2 — 1 — 2k 2 y 2 ) + diG 2 (y) — uj (31) 

V 2 (x) = ~ d 2 r 2 (2k 2 - 1 - 2k 2 y 2 ) + 4d 2 G 2 (y) - a 

Wj{pc) = jrdjdc(rx , k)sn(rx, k) ( 2 G(p) + yG'(y )) 

where j G { 1 , 2 }. Similar to the previous section, we only need to choose a function G(y) so that (2 G(y) + yG f (y )) is 
an even function of x in order to satisfy the VT symmetry criterion. Equations ( [29] ) and © give us a solution we 
call the linear oscillatory case. More details are included in Section [TV[ 


C. Other solutions 

Here, we introduce a possibility which is distinct from the previous ones as follows. We introduce F\(y) and F 2 (y) 
in the form 


F, = iAy p ( 1 - y 2 ) 1 / 2 , F 2 = By* (32) 

for p, q > 0 i.e., a non-polynomial form. By © and examining conditions 0 -® we find that we must have 
p(p — 1)(1 — k 2 ) = q(q — 1)(1 — k 2 ) = 0 and q < 2 p. That is, we require that either k = 1 with q < 2p, or 
p = g = l, orp = g = 0. We will focus on the former two cases. As for in 0 and condition ([9|, we find that 
G(y) must be in the form 


G(y)=Cy a (l-y 2 ) b 


(33) 


where a, b G IN. 


i. Cnoidal parameter k — 1 and y — sech(ra;) 


In the k = 1 case, we have the solutions 


cj = -dir 2 p 2 , 
cr = -d 2 r 2 q 2 . 


Vi{x) = -By q + d 1 r 2 y 2 (p + l)(p + 2) + d 1 G 2 (y) 

V 2 (x) = -T y 2 P - 9(1 _ y2) + d 2 V 2 y‘ 2 q(q + 1) + 4d 2 G 2 (y) 


(34) 


with q < 2 p and G(p) in the form of 
clear which choices of a;, cr will give Vq2 
made in (34). Also, we have by 0 


Since for this family of solutions the form of G is specified, it is immediately 
0 as x oo. In contrast to previous sections, those choices have been 


W\ = Grdisech(rx) tanh(rx) [y a 1 (1 — y 2 ) b (a + 2 p) — 2 p a+1 (l — y 2 ) h 1 (6 + 1 )) 

W 2 = 2Grd 2 tdmh(rx)(y a (l-y 2 ) b (a-h2q)-2by a+2 (l-y 2 ) b - 1 ). (35) 


Equations (32)-(|35|) give us a solution that bears a bright soliton coupled with a dark-in-bright soliton. The latter 
involves a pair of bright solitary waves coupled in a bound state anti-symmetric (i.e., they bear a phase difference 
of 7 r) configuration; another example of this form has been previously reported e.g. in m • More details on the 
propagation of this solution and its stability are included in Section m 


2. Powers p = q = 1 


In the case where p = q = 1, we obtain 
ijj = — dir 2 (bk 2 — 4), 
a = —d 2 r 2 (2k 2 — 1 ), 


Vi(x) = —By + dd\r 2 k 2 y 2 + diG 2 (y) 

V 2 {x) = -T y (i _ y 2 ) + 2 d 2 r 2 k 2 y 2 + 4 d 2 G 2 (y). 


( 36 ) 




Again here since G is known we have made choices of cj,ct reflected in (36) so that Vi^ 0 as x oo. With G(y) 
as in (33) we have 


Wi = rCdisn(rx, k)dn(rx, k) (y a 1 (1 — y 2 ) b (a + 2) — 2y a+1 (l — y 2 ) b 1 (b- 
W 2 = 2 rCd 2 dc(rx , k)sn(rx, k) [y a { 1 — y 2 ) b (a + 2) — 2by a+2 (l — ?/ 2 ) 6_1 ) . 


1 )) 


(37) 


We will refer to the solution in (|32|-(|33j) and ([36j)-([37j) as the non-polynomial oscillatory solution. In the special case 
of k = 1, this reverts to a waveform of the same type as the one examined above (namely, a bright solitary wave 
coupled to a dark-in-bright one). More details are provided on this solution in Section [TV[ 


IV. STABILITY AND DYNAMICS OF THE SOLUTIONS 


To study the stability of solutions we will first present the corresponding linear stability analysis framework. We 
begin by writing 


* z\„—iujz 


u= (U(x) + a(x)e Xz + 6(cc)*e A z )e 


\* „\* z \ 2 iujz 


v = (y{x) + c(x)e Az + d(cc)*e A z )e 


(38) 


(u 

+ L\ 

-V 

-U* 

0 \ 

(<>\ 


/«' 

u* 

-LJ- L\ 

0 

u 

6 1 

= A 

b 


-2 U 

0 

2c d K> L 2 

0 

c 

c 

V ’ 

0 

2 U* 

0 

—2c j — hi — L 2 ) 

LJ 


\d. 


where U(x), V(x) are the exact solutions of © found in Section [TTT| Substituting ( [38] ) into the system § we obtain 
in the first order set of equations 


(39) 


where the operators L \ : 2 are Li = did xx + Vi(x) + iWi(x) for i G {1,2}. If, for a given solution, the corresponding 
eigenvalue A has a positive real part then the solution is unstable as is readily seen in (38); otherwise the solution is 
stable. 

In the following subsections, we will apply the linear stability analysis and show the results of numerical propagation 
of the solutions we found in Section III according to ([!]) using a standard explicit 4th order Runge-Kutta code. We 
focus on the three solitonic solutions derived in Section |III[ of dark-dark, bright-bright and also bright coupled with 
the dark-in-bright waveforms, and also on the three oscillatory solutions derived in Section [ITT] in each of the quadratic, 
linear and non-polynomial cases considered. In each subsection, we start by specifying G(y) and other parameters 
as is necessary. We find that in all cases the solutions are unstable with increasing strength of instability as the 
amplitude parameters increase. Each example we consider has various regions of weak and strong instability as is 
discussed in the following subsections. 

Notice also that the equation 




! J {W\{x)\u(x, z)\ 2 + W 2 (x)\v(x, z)\ 2 ) dx 


(40) 


can be derived from the system (HT) where the combined power function P(z ) is defined as P(z) = 
J [\u{x , z)\ 2 + \v(x, z)| 2 ) dx. Equation (|40|) acts as a numerical check of all of the simulations performed in this 


section. 


A. Solitonic Solutions (k = 1) 

1. Dark-dark solitary wave 


For the solutions presented in Section [III A| in equations (|2d|)-(|2l|) we additionally make the choice here of G(y) = 

(41) 


KFi(y) with K G R. Also choosing a;, a so that » 0 as x 00 gives 

3 


U} = ACo + diK 2 , Vi = ^4Co(l - ^y 2 ) + 6dir 2 y 2 + diK 2 {l - -y 2 ) 2 - lo 

a = < L+4d 2 K 2 , V 2 = t}(l - ^y 2 ) + 6d 2 r 2 y 2 +4d 2 K 2 (l - ^y 2 ) 2 - a 


Wj = —9jrdjKtdinh.(rx)sech 2 (rx). 


( 42 ) 
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FIG. 1: This figure shows stability analysis and time propagation plots for the dark-dark soliton case. The parameters 
di = d ,2 = 0.1, r = 1 and K — 0.05 are fixed in every plot. The contour plot in the upper left depicts max(Re(A)) as a function 
of amplitude parameters Co and ACo of the F\, functions. The other seven plots correspond to the point (1,1) in the dim 
but non-zero region of the contour plot. At this point we have max(Re(A)) ~ 0.5201. In the left two plots of the bottom 
row, we show the magnitudes of the real and imaginary parts of the potential functions: \Vi(x)\, |V 2 | and \Wi(x)\, |W 2 (^)| 
respectively. In the top two plots of the center column, we show the magnitudes of the eigenvectors \a(x)\ (blue), \c(x)\ (green) 
and the eigenvalues A in the complex plane that correspond to the stationary solutions seen in the top right panel. The right 
column shows the magnitudes of the solution at t = 0 and at later times. Here we find that the unstable solution loses its dark 
soliton shape over time, with the destabilization manifesting across much of the x axis. 


We present the stability analysis of this family of solutions in Figure [l] We find that as the amplitudes Co, ACo of the 
F u F 2 functions, respectively, increase the solution becomes increasingly unstable. The panels of time propagation 
plots in Figure [l] show that over the dynamical evolution, these unstable dark soliton solutions will not maintain 
the dark soliton shape. Instead, the wide range of unstable eigenmodes in the system will induce a form of “lattice 
turbulence” whereby the end dynamical result will appear to bear no clear solitonic (or other) structure. 


2. Antidark-antidark solitary wave 


For the solutions presented in Section III A in equations (22)-(23) with the 
K G IR. Also choosing lj, a so that W ,2 0 as x —>> 00 gives 


sign, we take G(y) = KF\(y) with 


u =AC 0 J rd 1 K 2 , 
a =^+4 d 2 K 2 , 


Vi = AC 0 (V22y + Ay 2 ) - d l r 2 {V 1 22y - 6 y 2 ) + ckK 2 ( 1 + V22y + Ay 2 ) 2 - d x K 2 (43) 

r 2 = *=j(V22y + Ay 2 ) - d 2 r 2 (V22y - 6 y 2 ) + Ad 2 K 2 {l + V22y + Ay 2 ) 2 - Ad 2 K 2 
Wj = 3jrdjKtdinh(rx)sech(rx) (V‘22 + 8sech(rx) J . (44) 


We present the stability analysis of this family of solutions in Figure |2j We find that as the amplitudes Co,ALCo 
of the Fi,F 2 functions, respectively, increase the solution becomes increasingly unstable with a pattern similar to 
the previous example. However, here the eigenvectors are localized. The panels of time propagation plots in Figure 
[2] show that over the dynamical evolution, these soliton solutions will not maintain the soliton shape. Instead, the 
turbelence occurs near the center of the lattice, close to the solution’s peak. The instability is similar over time to 
what is observed in the dark-in-bright example below. 
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FIG. 2: This figure shows stability analysis and time propagation plots for the antidark-antidark soliton case. The placement 
of the figures follows the same pattern as that of Figure [l] The parameters di=c ?2 = 0.1, r = l and K — 0.05 are fixed in 
every plot. The point C — D = 0.5 from the upper left contour plot corresponds to max(Re(A)) ~ 1.5547, and the other seven 
plots show details about these amplitude values. Here we find that the unstable solution loses its soliton shape in a way that 
is similar to the dark-in-bright example below in Figure [5] 


3. Multiple dark solitary wave 


For the solutions presented in Section III A in equations (22)-(23) with the — sign, we take G(y) = KF\(y) with 
K G IR. Also choosing lj, a so that Vi ,2 0 as x —>> oo gives 


u) = AC 0 + d 1 K 2 , 
a = C ^.+M 2 K 2 , 


Vi = AC 0 {-V22y + 4 y 2 ) + d ir 2 {VZ2y + 6 y 2 ) + (hK 2 (l - V?2y + 4 y 2 ) 2 - d x K 2 (45) 
V 2 = ^j(-V22y + 4 y 2 ) + d 2 r 2 {V^y + 6 y 2 ) + M 2 K 2 {1 - V22y + Ay 2 ) 2 - Ad 2 K 2 
Wj = 3jrdjiFtanh(rx)sech(rx) (~V 22 + 8 sech (rx) \ . (46) 


We present the stability analysis of this family of solutions in Figure [3| We find that as the amplitudes Co, ACo of 
the Fi,F 2 functions, respectively, increase the solution becomes increasingly unstable with a pattern similar to the 
previous example. The panels of time propagation plots in Figure [3] show that over the dynamical evolution, these 
soliton solutions will not maintain the dark soliton shape. Instead, the turbelence occurs across the x axis. The 
instability is similar over time to what is observed in the first dark-dark example above. 


4- Bright-Bright solitary wave 


For the solutions presented in IIIB in equations (29) and (30) we simply choose G(y) = Ky with K G IR and we 
take cj, cr to be such that Vi ,2 —> 0 as x oo. This gives the solution 


UJ 

a 


—dir 2 pi 2 , Vi = ±Dy p2 + d\K 2 y 2 + dir 2 pi(pi + 1 )y 2 , W\ = rd\K{2pi + l)sech(rx)tanh(rx) (47) 

n 2 

-d 2 r 2 p 2 2 , V 2 = ±-—y 2pi ~ P2 + 4 d 2 K 2 y 2 + d 2 r 2 p 2 (p 2 + 1 )y 2 , W 2 = 2rd 2 K(2p 2 + l)sech(ra)tanh(ra). 


In Figure [4] we show the stability analysis for selected parameters. Similar to the dark-dark case, the solutions do 
not maintain their shape as time progresses and the strength of instability increases as both the multipliers C and 
D in ( [29] ) increase. Here, it is clear that the instability results in the breaking of the parity symmetry, leading to a 
symmetry-breaking pattern. 
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FIG. 3: This figure shows stability analysis and time propagation plots for the multiple dark soliton case. The placement of 
the figures follows the same pattern as that of Figure^ The parameters d\ — cfe = 0.1, r = 1 and K — 0.05 are fixed in every 
plot. The point C — D — 0.5 from the upper left contour plot corresponds to max(Re(A)) ~ 0.5277, and the other seven plots 
show details about these amplitude values. Here we find that the unstable solution loses its dark soliton shape similar to the 
first dark-dark example. 



FIG. 4: This figure shows stability analysis and time propagation plots for the bright-bright soliton case with p = q = 1. The 
placement of the figures follows the same pattern as that of Figure [l] Here the parameters common to all plots are d\ — cfo = 1, 
r — 1 and K — 0.1. The point C — 0.5, D — 3 from the upper left contour plot corresponds to max(Re(A)) ~ 1.7055 and more 
details regarding these specific amplitude parameters are shown in the other seven plots. Here we find that as time progresses 
the solution loses its soliton shape with turbulent, symmetry-breaking behaviour occurring first nearby the central peak and 
then leaking outward across the x axis. 
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FIG. 5: This figure shows stability analysis and time propagation plots for the bright and dark-in-bright soliton case. The 
placement of the figures follows the same pattern as that of Figure [l] Here the parameters common to all plots are d\ — cfo = 0.1, 
r = 0.2 and K = 0.05. The point C — 0.5, D — 3 from the upper left contour plot corresponds to max(Re(A)) ~ 0.2778 and 
more details regarding these specific amplitude parameters are shown in the other seven plots. Here we find that as time 
progresses the solution loses its soliton shape completely. 


5. Bright and Dark-in-Bright solitary Wave 


For the solutions presented in Section III C in equations 
lj = —dir 2 p 2 , 
a = -d 2 r 2 g 2 , 


35| we choose G(y) = Ky( 1 — y 2 ) and obtain 


,. 2\2 


Vi = -By q + d\r 2 y 2 (p + l)(p + 2) + d x K 2 y 2 (\ - y 2 ) 

V 2 = ~^-y 2p ~ g ( 1 - y 2 ) + d 2 r 2 y 2 q{q + 1) + 4d 2 K 2 y 2 (l - y 2 ) 2 

Wi = Krdisech(rx) tanh(rx)(l + 2p — (5 + 2 p)y 2 ) 

W 2 = 2Krd2sech(rx) tanh(rx)(l + 2q — (3 + 2 q)y 2 ) 


by choosing as usual. The stability analysis is presented in Figure [5] Similar to the quadratic case, the strength 
of instability increases as A increases and as B increases. The propagation plots in Figure [5] show that the peak 
destabilizes and the amplitude spreads out over the x axis while maintaining some comparative concentration at the 
center of the axis. Furthermore, a symmetry-breaking feature appears once again to be amplifed and be distinctly 
observable at the end of the simulation’s reporting horizon. 


B. Oscillatory Solutions (k = 0) 

1. Quadratic oscillatory solution 


For the solutions in Section III A in equations 
case, obtaining the following solutions 

ijj = ACq + d\K 2 + 4dir 2 , 

a = 9± + M 2 K 2 + M 2r 2 , 


(27)-(28) we take G(y) 


KFi(y) similar to the dark-dark soliton 


Vi = -2 AC 0 y 2 + diif 2 (l - 2y 2 ) 2 - diK 2 (48) 

V 2 = -2^V + 4d 2 K 2 (l - 2y 2 ) 2 - 4d 2 K 2 
Wj = —12 jrdjK cos (rx) sin (rx) 


for j G {1,2}. The stability graph in Figure [ 6 ] has similar features to the one in Figure [I] showing that the changes 
in stability strength of the system across the amplitudes Co, ACq grid are similar despite very different k values. 
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FIG. 6: This figure shows stability analysis and time propagation plots for the quadratic oscillatory case. The placement of the 
figures follows the same pattern as that of Figure [l] Here the parameters common to all plots are d± = d 2 = 0.1, r = 0.5 and 
K = 0.05. The point C — 0.5, D — 3 from the upper left contour plot corresponds to max(Re(A)) ~ 2.3602 and more details 
regarding these specific amplitude parameters are shown in the other seven plots. Here we find that as time progresses the 
instability of the solution appears across the x axis at near regular intervals, close to the periodicity of the stationary solution. 
However, the solution over time does not remain truly periodic. 


In the current oscillatory function case, we observe that the waves will not maintain their original shapes, with the 
most apparent distortions located at near-periodic points along the x axis. These distortions will not only break the 
periodicity of the structure but they will also lead (within some lattice periods of the solution) into the turbulent 
dynamical evolution discussed previously. 


2. Linear oscillatory solution 


For the solutions in Section IIIB in equations (29) and (31) we take G(y) = Ky and obtain the solution 


c u = dir 2 , V\ = ±Dy + d\K 2 y 2 , W\ = 3rd\K sin(rx) 

C 2 

a = d 2 r 2 , V 2 = ±—y + ^d 2 K 2 y 2 , W 2 = 6rd 2 K sin (rx). 


(49) 


The stability graph in this linear function case for k = 0 is presented in Figure [7| Here we see that the pattern of 
the strength of the instability is more similar to that of the quadratic functions than it is to the linear functions case 
with k = 1. We can see that the strength of instability increases as C, D increase. The propagation plots in Figure [7] 
show that over time the wave loses its original shape at points across the x axis. Here the instability is induced by 
eigenvectors which also extend across the x axis and which lead to a breakup of the periodicity of the original pattern. 


3. Other oscillatory solution 


For the solutions in Section IIIC in equations (32)-(33) and (36)-(37) we choose G(y) = Ky{l — y 2 ) 
and obtain 


= a = 1 , 


p = q 


uj =4dir 2 , V 1 = -By + d 1 K 2 y 2 (l-y 2 ) 2 (50) 

a =d 2 r 2 , V 2 = -^-y(l - y 2 ) +M 2 K 2 y 2 (l - y 2 ) 2 

W\ = Krd\ sin {rx) (3 — 7 cos 2 (rx)) 

W 2 = 2Krd 2 sin(rx) (3 — 5 cos 2 (rx)) . 
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FIG. 7: This figure shows stability analysis and time propagation plots for the linear oscillatory case. The placement of the 
figures follows the same pattern as that of Figure [l] Here the parameters common to all plots are d\ — d<i = 0.35, r = 0.25 and 
K = 0.1. The point C = 0.5, D = 3 from the upper left contour plot corresponds to max(Re(A)) ~ 2.6630 and more details 
regarding these specific amplitude parameters are shown in the other seven plots. Similar to the quadratic oscillatory case, 
here we find that as time progresses the solution does not remain truly periodic. 
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FIG. 8: This figure shows stability analysis and time propagation plots for the other oscillatory case. The placement of the 
figures follows the same pattern as that of Figure^ Here the parameters common to all plots are d\ — cfo = 0.1, r = 0.25 and 
K — 0.05. The point C — 0.5, D — 3 from the upper left contour plot corresponds to max(Re(A)) ~ 2.8819 and more details 
regarding these specific amplitude parameters are shown in the other seven plots. Here we find that as time progresses the 
instability of the solution appears across the x axis at near regular intervals, close to the periodicity of the stationary solution. 
The instability is similar to the other oscillatory cases and again leads to a periodicity breakup. 


The stability graph on the left of Figure [8] shows that the strength of instability increases as B increases, yet it appears 
to be roughly independent of A. The propagation panels show the distortion of the original solution occurring over 
time at points across the x axis corresponding to an eigenvector that is also spread across x. 





































15 


V. CONCLUSION 

In the present work, we have explored both solitary and more broadly periodic (including cnoidal and even their 
trigonometric limit of k = 0) solutions of the VT -symmetric problem with quadratic nonlinearity. A reverse engineer¬ 
ing approach was adopted herein attempting to identify even real potentials and odd imaginary ones that would be 
compatible with specific cnoidal solutions (and their hyperbolic limits in the case of k = 1, as well as their trigono¬ 
metric ones in the case of k = 0). It was shown that necessitating the existence of such solutions generally leads 
to a number of plausible requirements (for the absence of singularities) that can, in turn, be used to identify wide 
parametric families of potentials with the desired solutions. Relevant waveforms included, but were arguably not 
limited to dark-dark or bright-bright solitary waves and more exotic generalizations thereof such as the bright wave 
coupled to a dark-in-bright structure. Oscillatory variants of such hyperbolic limit solutions were identified as well. 

Naturally, numerous directions of future research arise from the present considerations. Offering a systematic similar 
approach could be of interest also in the case of other nonlinearities. From a stability perspective, it would appear 
interesting to identify case examples with stable isolated parameter values or, more promisingly, wide parameter 
ranges, as the solutions considered here seemed to be largely unstable (with bands of unstable modes) resulting in 
turbulent dynamics in many of our dynamical examples. Finally, exploring two-dimensional generalizations of the 
relevant PT-symmetric systems is of particular interest in its own right both at the level of discrete systems (see e.g. 
the plaquette considerations of M) and at that of continuum ones (see e.g. 0 ); see also the recent work of [52]. Such 
studies are currently in progress and will be reported in future publications. 
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